Regional association plotting of 11 loci associated with CAC.

A ‘druggable-MI-targets’ project

Setup

We will clean the environment, setup the locations, define colors, and create a datestamp.

Clean the environment.

Set locations and working directories…

Create a new analysis directory...
[1] FALSE
[1] FALSE
[1] FALSE
[1] FALSE
[1] FALSE
[1] FALSE
[1] "/Users/swvanderlaan/PLINK/analyses/consortia/CHARGE_1000G_CAC"
 [1] "1. CHARGE_1000G_CAC.nb.html"             "1. CHARGE_1000G_CAC.Rmd"                 "2. bulkRNAseq.nb.html"                   "2. bulkRNAseq.Rmd"                       "20220125.CAC.Parsing_GWASSumStats.RData"
 [6] "3. scRNAseq.nb.html"                     "3. scRNAseq.Rmd"                         "4. Parsing_GWASSumStats.nb.html"         "4. Parsing_GWASSumStats.Rmd"             "5. RegionalAssociationPlots.nb.html"    
[11] "5. RegionalAssociationPlots.Rmd"         "CAC"                                     "CHARGE_1000G_CAC.Rproj"                  "CredibleSets"                            "images"                                 
[16] "LICENSE"                                 "RACER"                                   "README.html"                             "README.md"                               "README.orig.md"                         
[21] "renv"                                    "renv.lock"                               "scripts"                                 "SNP"                                     "targets"                                

… a package-installation function …

source(paste0(PROJECT_loc, "/scripts/functions.R"))

… and load those packages.

install.packages.auto("readr")

Loading required package: readr

install.packages.auto("optparse")

Loading required package: optparse

install.packages.auto("tools")

Loading required package: tools

install.packages.auto("dplyr")

Loading required package: dplyr

Attaching package: 'dplyr'

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

install.packages.auto("tidyr")

Loading required package: tidyr

install.packages.auto("naniar")

Loading required package: naniar

# To get 'data.table' with 'fwrite' to be able to directly write gzipped-files
# Ref: https://stackoverflow.com/questions/42788401/is-possible-to-use-fwrite-from-data-table-with-gzfile
# install.packages("data.table", repos = "https://Rdatatable.gitlab.io/data.table")
library(data.table)

data.table 1.14.2 using 1 threads (see ?getDTthreads).  Latest news: r-datatable.com
**********
This installation of data.table has not detected OpenMP support. It should still work but in single-threaded mode.
This is a Mac. Please read https://mac.r-project.org/openmp/. Please engage with Apple and ask them for support. Check r-datatable.com for updates, and our Mac instructions here: https://github.com/Rdatatable/data.table/wiki/Installation. After several years of many reports of installation problems on Mac, it's time to gingerly point out that there have been no similar problems on Windows or Linux.
**********

Attaching package: 'data.table'

The following objects are masked from 'package:dplyr':

    between, first, last

install.packages.auto("tidyverse")

Loading required package: tidyverse
── Attaching packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.1 ──
✓ ggplot2 3.3.5     ✓ stringr 1.4.0
✓ tibble  3.1.6     ✓ forcats 0.5.1
✓ purrr   0.3.4     
── Conflicts ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
x data.table::between() masks dplyr::between()
x dplyr::filter()       masks stats::filter()
x data.table::first()   masks dplyr::first()
x dplyr::lag()          masks stats::lag()
x data.table::last()    masks dplyr::last()
x purrr::transpose()    masks data.table::transpose()

install.packages.auto("knitr")

Loading required package: knitr

install.packages.auto("DT")

Loading required package: DT

install.packages.auto("eeptools")

Loading required package: eeptools
Welcome to eeptools for R version 1.2.0!
Developed by Jared E. Knowles 2012-2018
for the Wisconsin Department of Public Instruction
Distributed without warranty.

install.packages.auto("haven")

Loading required package: haven

install.packages.auto("tableone")

Loading required package: tableone

install.packages.auto("BlandAltmanLeh")

Loading required package: BlandAltmanLeh

# Install the devtools package from Hadley Wickham
install.packages.auto('devtools')

Loading required package: devtools
Loading required package: usethis

library(devtools) 

# for plotting
install.packages.auto("pheatmap")

Loading required package: pheatmap

install.packages.auto("forestplot")

Loading required package: forestplot
Loading required package: grid
Loading required package: magrittr

Attaching package: 'magrittr'

The following object is masked from 'package:purrr':

    set_names

The following object is masked from 'package:tidyr':

    extract

Loading required package: checkmate

install.packages.auto("ggplot2")
install.packages.auto("ggpubr")

Loading required package: ggpubr

install.packages.auto("ggrepel")

Loading required package: ggrepel

install.packages.auto("UpSetR")

Loading required package: UpSetR

devtools::install_github("thomasp85/patchwork")

Using github PAT from envvar GITHUB_PAT
Skipping install of 'patchwork' from a github remote, the SHA1 (79223d30) has not changed since last install.
  Use `force = TRUE` to force installation

# For regional association plots
install_github("oliviasabik/RACER") 

Using github PAT from envvar GITHUB_PAT
Skipping install of 'RACER' from a github remote, the SHA1 (1394c9d4) has not changed since last install.
  Use `force = TRUE` to force installation

# Install ggrepel package if needed
unloadNamespace('Seurat')
library(ggrepel)

We will create a datestamp and define the Utrecht Science Park Colour Scheme.

Today = format(as.Date(as.POSIXlt(Sys.time())), "%Y%m%d")
Today.Report = format(as.Date(as.POSIXlt(Sys.time())), "%A, %B %d, %Y")

### UtrechtScienceParkColoursScheme
###
### WebsitetoconvertHEXtoRGB:http://hex.colorrrs.com.
### Forsomefunctionsyoushoulddividethesenumbersby255.
###
### No. Color                 HEX   (RGB)                                     CHR         MAF/INFO
###---------------------------------------------------------------------------------------
### 1     yellow                #FBB820 (251,184,32)                      =>    1       or 1.0>INFO
### 2     gold                #F59D10 (245,157,16)                    =>    2       
### 3     salmon                #E55738 (229,87,56)                   =>    3       or 0.05<MAF<0.2 or 0.4<INFO<0.6
### 4     darkpink          #DB003F ((219,0,63)                   =>    4       
### 5     lightpink         #E35493 (227,84,147)                      =>    5       or 0.8<INFO<1.0
### 6     pink                #D5267B (213,38,123)                    =>    6       
### 7     hardpink          #CC0071 (204,0,113)                   =>    7       
### 8     lightpurple       #A8448A (168,68,138)                      =>    8       
### 9     purple                #9A3480 (154,52,128)                      =>    9       
### 10  lavendel            #8D5B9A (141,91,154)                      =>    10      
### 11  bluepurple        #705296 (112,82,150)                    =>    11      
### 12  purpleblue        #686AA9 (104,106,169)               =>    12      
### 13  lightpurpleblue #6173AD (97,115,173/101,120,180)    =>  13      
### 14  seablue             #4C81BF (76,129,191)                      =>    14      
### 15  skyblue             #2F8BC9 (47,139,201)                      =>    15      
### 16  azurblue            #1290D9 (18,144,217)                      =>    16      or 0.01<MAF<0.05 or 0.2<INFO<0.4
### 17  lightazurblue     #1396D8 (19,150,216)                    =>    17      
### 18  greenblue           #15A6C1 (21,166,193)                      =>    18      
### 19  seaweedgreen      #5EB17F (94,177,127)                    =>    19      
### 20  yellowgreen       #86B833 (134,184,51)                    =>    20      
### 21  lightmossgreen  #C5D220 (197,210,32)                      =>    21      
### 22  mossgreen           #9FC228 (159,194,40)                      =>    22      or MAF>0.20 or 0.6<INFO<0.8
### 23  lightgreen      #78B113 (120,177,19)                      =>    23/X
### 24  green                 #49A01D (73,160,29)                     =>    24/Y
### 25  grey                  #595A5C (89,90,92)                        =>  25/XY   or MAF<0.01 or 0.0<INFO<0.2
### 26  lightgrey           #A2A3A4 (162,163,164)                 =>    26/MT
###
### ADDITIONAL COLORS
### 27  midgrey         #D7D8D7
### 28  verylightgrey   #ECECEC"
### 29  white           #FFFFFF
### 30  black           #000000
###----------------------------------------------------------------------------------------------

uithof_color = c("#FBB820","#F59D10","#E55738","#DB003F","#E35493","#D5267B",
                 "#CC0071","#A8448A","#9A3480","#8D5B9A","#705296","#686AA9",
                 "#6173AD","#4C81BF","#2F8BC9","#1290D9","#1396D8","#15A6C1",
                 "#5EB17F","#86B833","#C5D220","#9FC228","#78B113","#49A01D",
                 "#595A5C","#A2A3A4", "#D7D8D7", "#ECECEC", "#FFFFFF", "#000000")

uithof_color_legend = c("#FBB820", "#F59D10", "#E55738", "#DB003F", "#E35493",
                        "#D5267B", "#CC0071", "#A8448A", "#9A3480", "#8D5B9A",
                        "#705296", "#686AA9", "#6173AD", "#4C81BF", "#2F8BC9",
                        "#1290D9", "#1396D8", "#15A6C1", "#5EB17F", "#86B833",
                        "#C5D220", "#9FC228", "#78B113", "#49A01D", "#595A5C",
                        "#A2A3A4", "#D7D8D7", "#ECECEC", "#FFFFFF", "#000000")
### ----------------------------------------------------------------------------

Introduction

We will parse the data to create regional association plots for each of the 11 loci.

Load data

We need to load the data first.

gwas_sumstats_racer <- readRDS(file = paste0(OUT_loc, "/gwas_sumstats_racer.rds"))

Regional association plotting

Top 11 loci

We are interested in 11 top loci.

library(openxlsx)
variant_list <- read.xlsx(paste0(TARGET_loc, "/Variants.xlsx"), sheet = "TopLoci")

DT::datatable(variant_list)

NA

Let’s do some plotting.

library(RACER)
# Make directory for plots
ifelse(!dir.exists(file.path(PROJECT_loc, "/RACER")), 
       dir.create(file.path(PROJECT_loc, "/RACER")), 
       FALSE)

[1] FALSE

RACER_loc = paste0(PROJECT_loc,"/RACER")

variants_of_interest <- c(variant_list$rsID)

variants_of_interest_fewgenes <- c("rs9349379", "rs3844006", "rs2854746", "rs4977575", "rs9633535", "rs11063120", "rs9515203", "rs7182103")

for(VARIANT in variants_of_interest){
  cat(paste0("Getting data for ", VARIANT,".\n"))

  tempCHR <- subset(variant_list, rsID == VARIANT)[,5]
  tempSTART <- subset(variant_list, rsID == VARIANT)[,17]
  tempEND <- subset(variant_list, rsID == VARIANT)[,18]
  tempVARIANTnr <- subset(variant_list, rsID == VARIANT)[,1]

  cat("\nSubset required data.\n")
  temp <- subset(gwas_sumstats_racer, Chr == tempCHR & (Position >= tempSTART & Position <= tempEND))
  
  cat("\nFormatting association data.\n")
  temp_f = RACER::formatRACER(assoc_data = temp, chr_col = 3, pos_col = 4, p_col = 5)

  cat("\nGetting LD data.\n")
  temp_f_ld = RACER::ldRACER(assoc_data = temp_f, rs_col = 2, pops = "EUR", lead_snp = VARIANT)
  
  cat(paste0("\nPlotting region surrounding ", VARIANT," on ",tempCHR,":",tempSTART,"-",tempEND,".\n"))
  # source(paste0(PROJECT_loc, "/scripts/functions.R"))
  p1 <- singlePlotRACER2(assoc_data = temp_f_ld, 
                               chr = tempCHR, build = "hg19", 
                               plotby = "snp", snp_plot = VARIANT,
                               label_lead = TRUE, gene_track_h = 2, gene_name_s = 1.75)
  
  print(p1)
  cat(paste0("Saving image for ", VARIANT,".\n"))
  # ggsave(filename = paste0(RACER_loc, "/", tempVARIANTnr, ".", Today, ".",VARIANT,".regional_assoc.png"), plot = last_plot())
  # ggsave(filename = paste0(RACER_loc, "/", tempVARIANTnr, ".", Today, ".",VARIANT,".regional_assoc.pdf"), plot = last_plot())
  ggsave(filename = paste0(RACER_loc, "/", tempVARIANTnr, ".", Today, ".",VARIANT,".regional_assoc.eps"), plot = last_plot())
  
  rm(temp, p1,
     temp_f, temp_f_ld,
     tempCHR, tempSTART, tempEND,
     VARIANT, tempVARIANTnr)
  
}

Getting data for rs9349379.

Subset required data.

Formatting association data.

Formating association data...
Processing -log10(p-values)...
Preparing association data...

Getting LD data.

All inputs are go!
Reading in association data...
Populations selected.
Calculating LD using rs9349379...

[1] "https://ldlink.nci.nih.gov/LDlinkRest/ldproxy?var=rs9349379&pop=EUR&r2_d=r2&token=c0f613f149ab"

 [100%] Downloaded 70967 bytes...
Merging input association data with LD...

Plotting region surrounding rs9349379 on 6:12403957-13403957.

Plotting by...
snp rs9349379
Reading in association data
Generating Plot

Saving image for rs9349379.

Saving 7.29 x 4.51 in image

Getting data for rs3844006.

Subset required data.

Formatting association data.

Formating association data...
Processing -log10(p-values)...
Preparing association data...

Getting LD data.

All inputs are go!
Reading in association data...
Populations selected.
Calculating LD using rs3844006...

[1] "https://ldlink.nci.nih.gov/LDlinkRest/ldproxy?var=rs3844006&pop=EUR&r2_d=r2&token=c0f613f149ab"

 [100%] Downloaded 31016 bytes...
Merging input association data with LD...

Plotting region surrounding rs3844006 on 6:131595002-132595002.

Plotting by...
snp rs3844006
Reading in association data
Generating Plot

Saving image for rs3844006.

Saving 7.29 x 4.51 in image

Getting data for rs2854746.

Subset required data.

Formatting association data.

Formating association data...
Processing -log10(p-values)...
Preparing association data...

Getting LD data.

All inputs are go!
Reading in association data...
Populations selected.
Calculating LD using rs2854746...

[1] "https://ldlink.nci.nih.gov/LDlinkRest/ldproxy?var=rs2854746&pop=EUR&r2_d=r2&token=c0f613f149ab"

 [44%] Downloaded 48575 bytes...
 [100%] Downloaded 109298 bytes...
Merging input association data with LD...

Plotting region surrounding rs2854746 on 7:45460645-46460645.

Plotting by...
snp rs2854746
Reading in association data
Generating Plot

Saving image for rs2854746.

Saving 7.29 x 4.51 in image

Getting data for rs4977575.

Subset required data.

Formatting association data.

Formating association data...
Processing -log10(p-values)...
Preparing association data...

Getting LD data.

All inputs are go!
Reading in association data...
Populations selected.
Calculating LD using rs4977575...

[1] "https://ldlink.nci.nih.gov/LDlinkRest/ldproxy?var=rs4977575&pop=EUR&r2_d=r2&token=c0f613f149ab"

 [100%] Downloaded 26000 bytes...
Merging input association data with LD...

Plotting region surrounding rs4977575 on 9:21624744-22624744.

Plotting by...
snp rs4977575
Reading in association data
Generating Plot

Saving image for rs4977575.

Saving 7.29 x 4.51 in image

Getting data for rs10899970.

Subset required data.

Formatting association data.

Formating association data...
Processing -log10(p-values)...
Preparing association data...

Getting LD data.

All inputs are go!
Reading in association data...
Populations selected.
Calculating LD using rs10899970...

[1] "https://ldlink.nci.nih.gov/LDlinkRest/ldproxy?var=rs10899970&pop=EUR&r2_d=r2&token=c0f613f149ab"

 [30%] Downloaded 48575 bytes...
 [71%] Downloaded 114111 bytes...
 [91%] Downloaded 146879 bytes...
 [100%] Downloaded 160476 bytes...
Merging input association data with LD...

Plotting region surrounding rs10899970 on 10:44015716-45334720.

Plotting by...
snp rs10899970
Reading in association data
Generating Plot

Saving image for rs10899970.

Saving 7.29 x 4.51 in image

Getting data for rs9633535.

Subset required data.

Formatting association data.

Formating association data...
Processing -log10(p-values)...
Preparing association data...

Getting LD data.

All inputs are go!
Reading in association data...
Populations selected.
Calculating LD using rs9633535...

[1] "https://ldlink.nci.nih.gov/LDlinkRest/ldproxy?var=rs9633535&pop=EUR&r2_d=r2&token=c0f613f149ab"

 [100%] Downloaded 49363 bytes...
Merging input association data with LD...

Plotting region surrounding rs9633535 on 10:63336088-64336088.

Plotting by...
snp rs9633535
Reading in association data
Generating Plot

Saving image for rs9633535.

Saving 7.29 x 4.51 in image

Getting data for rs10762577.

Subset required data.

Formatting association data.

Formating association data...
Processing -log10(p-values)...
Preparing association data...

Getting LD data.

All inputs are go!
Reading in association data...
Populations selected.
Calculating LD using rs10762577...

[1] "https://ldlink.nci.nih.gov/LDlinkRest/ldproxy?var=rs10762577&pop=EUR&r2_d=r2&token=c0f613f149ab"

 [100%] Downloaded 92997 bytes...
Merging input association data with LD...

Plotting region surrounding rs10762577 on 10:75417431-76417431.

Plotting by...
snp rs10762577
Reading in association data
Generating Plot

Saving image for rs10762577.

Saving 7.29 x 4.51 in image

Getting data for rs11063120.

Subset required data.

Formatting association data.

Formating association data...
Processing -log10(p-values)...
Preparing association data...

Getting LD data.

All inputs are go!
Reading in association data...
Populations selected.
Calculating LD using rs11063120...

[1] "https://ldlink.nci.nih.gov/LDlinkRest/ldproxy?var=rs11063120&pop=EUR&r2_d=r2&token=c0f613f149ab"

 [100%] Downloaded 41883 bytes...
Merging input association data with LD...

Plotting region surrounding rs11063120 on 12:3986618-4986618.

Plotting by...
snp rs11063120
Reading in association data
Generating Plot

Saving image for rs11063120.

Saving 7.29 x 4.51 in image

Getting data for rs9515203.

Subset required data.

Formatting association data.

Formating association data...
Processing -log10(p-values)...
Preparing association data...

Getting LD data.

All inputs are go!
Reading in association data...
Populations selected.
Calculating LD using rs9515203...

[1] "https://ldlink.nci.nih.gov/LDlinkRest/ldproxy?var=rs9515203&pop=EUR&r2_d=r2&token=c0f613f149ab"

 [100%] Downloaded 30139 bytes...
Merging input association data with LD...

Plotting region surrounding rs9515203 on 13:110549623-111549623.

Plotting by...
snp rs9515203
Reading in association data
Generating Plot

Saving image for rs9515203.

Saving 7.29 x 4.51 in image

Getting data for rs7182103.

Subset required data.

Formatting association data.

Formating association data...
Processing -log10(p-values)...
Preparing association data...

Getting LD data.

All inputs are go!
Reading in association data...
Populations selected.
Calculating LD using rs7182103...

[1] "https://ldlink.nci.nih.gov/LDlinkRest/ldproxy?var=rs7182103&pop=EUR&r2_d=r2&token=c0f613f149ab"

 [41%] Downloaded 48575 bytes...
 [82%] Downloaded 97745 bytes...
 [100%] Downloaded 118299 bytes...
Merging input association data with LD...

Plotting region surrounding rs7182103 on 15:78623946-79623946.

Plotting by...
snp rs7182103
Reading in association data
Generating Plot

Saving image for rs7182103.

Saving 7.29 x 4.51 in image

Getting data for rs7412.

Subset required data.

Formatting association data.

Formating association data...
Processing -log10(p-values)...
Preparing association data...

Getting LD data.

All inputs are go!
Reading in association data...
Populations selected.
Calculating LD using rs7412...

[1] "https://ldlink.nci.nih.gov/LDlinkRest/ldproxy?var=rs7412&pop=EUR&r2_d=r2&token=c0f613f149ab"

 [100%] Downloaded 38865 bytes...
Merging input association data with LD...

Plotting region surrounding rs7412 on 19:44912079-45912079.

Plotting by...
snp rs7412
Reading in association data
Generating Plot

Saving image for rs7412.

Saving 7.29 x 4.51 in image

variants_of_interest_manygenes <- c("rs7412", "rs10762577")
source(paste0(PROJECT_loc, "/scripts/functions.R"))

for(VARIANT in variants_of_interest_manygenes){
  cat(paste0("Getting data for ", VARIANT,".\n"))

  tempCHR <- subset(variant_list, rsID == VARIANT)[,5]
  tempSTART <- subset(variant_list, rsID == VARIANT)[,17]
  tempEND <- subset(variant_list, rsID == VARIANT)[,18]
  tempVARIANTnr <- subset(variant_list, rsID == VARIANT)[,1]

  cat("\nSubset required data.\n")
  temp <- subset(gwas_sumstats_racer, Chr == tempCHR & (Position >= tempSTART & Position <= tempEND))
  
  cat("\nFormatting association data.\n")
  temp_f = RACER::formatRACER(assoc_data = temp, chr_col = 3, pos_col = 4, p_col = 5)

  cat("\nGetting LD data.\n")
  temp_f_ld = RACER::ldRACER(assoc_data = temp_f, rs_col = 2, pops = "EUR", lead_snp = VARIANT)
  
  cat(paste0("\nPlotting region surrounding ", VARIANT," on ",tempCHR,":",tempSTART,"-",tempEND,".\n"))
  p1 <- singlePlotRACER2(assoc_data = temp_f_ld, 
                               chr = tempCHR, build = "hg19", 
                               plotby = "snp", snp_plot = VARIANT,
                               label_lead = TRUE, gene_track_h = 0.75, gene_name_s = 1.75)
  
  print(p1)
  cat(paste0("Saving image for ", VARIANT,".\n"))
  # ggsave(filename = paste0(RACER_loc, "/", tempVARIANTnr, ".", Today, ".",VARIANT,".regional_assoc.png"), plot = last_plot())
  # ggsave(filename = paste0(RACER_loc, "/", tempVARIANTnr, ".", Today, ".",VARIANT,".regional_assoc.pdf"), plot = last_plot())
  ggsave(filename = paste0(RACER_loc, "/", tempVARIANTnr, ".", Today, ".",VARIANT,".regional_assoc.eps"), plot = last_plot())
  
  rm(temp, p1,
     temp_f, temp_f_ld,
     tempCHR, tempSTART, tempEND,
     VARIANT, tempVARIANTnr)
  
}

Getting data for rs7412.

Subset required data.

Formatting association data.

Formating association data...
Processing -log10(p-values)...
Preparing association data...

Getting LD data.

All inputs are go!
Reading in association data...
Populations selected.
Calculating LD using rs7412...

[1] "https://ldlink.nci.nih.gov/LDlinkRest/ldproxy?var=rs7412&pop=EUR&r2_d=r2&token=c0f613f149ab"

 [100%] Downloaded 38865 bytes...
Merging input association data with LD...

Plotting region surrounding rs7412 on 19:44912079-45912079.

Plotting by...
snp rs7412
Reading in association data
Generating Plot

Saving image for rs7412.

Saving 7.29 x 4.51 in image

Getting data for rs10762577.

Subset required data.

Formatting association data.

Formating association data...
Processing -log10(p-values)...
Preparing association data...

Getting LD data.

All inputs are go!
Reading in association data...
Populations selected.
Calculating LD using rs10762577...

[1] "https://ldlink.nci.nih.gov/LDlinkRest/ldproxy?var=rs10762577&pop=EUR&r2_d=r2&token=c0f613f149ab"

 [100%] Downloaded 92997 bytes...
Merging input association data with LD...

Plotting region surrounding rs10762577 on 10:75417431-76417431.

Plotting by...
snp rs10762577
Reading in association data
Generating Plot

Saving image for rs10762577.

Saving 7.29 x 4.51 in image

variants_of_interest_cxcl12 <- c("rs10899970")
source(paste0(PROJECT_loc, "/scripts/functions.R"))

for(VARIANT in variants_of_interest_cxcl12){
  cat(paste0("Getting data for ", VARIANT,".\n"))

  tempCHR <- subset(variant_list, rsID == VARIANT)[,5]
  tempSTART <- subset(variant_list, rsID == VARIANT)[,17]
  tempEND <- subset(variant_list, rsID == VARIANT)[,18]
  tempVARIANTnr <- subset(variant_list, rsID == VARIANT)[,1]

  cat("\nSubset required data.\n")
  temp <- subset(gwas_sumstats_racer, Chr == tempCHR & (Position >= tempSTART & Position <= tempEND))
  
  cat("\nFormatting association data.\n")
  temp_f = RACER::formatRACER(assoc_data = temp, chr_col = 3, pos_col = 4, p_col = 5)

  cat("\nGetting LD data.\n")
  temp_f_ld = RACER::ldRACER(assoc_data = temp_f, rs_col = 2, pops = "EUR", lead_snp = VARIANT)
  
  cat(paste0("\nPlotting region surrounding ", VARIANT," on ",tempCHR,":",tempSTART,"-",tempEND,".\n"))
  p1 <- singlePlotRACER2(assoc_data = temp_f_ld, 
                               chr = tempCHR, build = "hg19", set = "all",
                               plotby = "snp", snp_plot = VARIANT,
                               label_lead = TRUE, gene_track_h = 0.75, gene_name_s = 1.75)
  
  print(p1)
  cat(paste0("Saving image for ", VARIANT,".\n"))
  # ggsave(filename = paste0(RACER_loc, "/", tempVARIANTnr, ".", Today, ".",VARIANT,".regional_assoc.png"), plot = last_plot())
  # ggsave(filename = paste0(RACER_loc, "/", tempVARIANTnr, ".", Today, ".",VARIANT,".regional_assoc.pdf"), plot = last_plot())
  ggsave(filename = paste0(RACER_loc, "/", tempVARIANTnr, ".", Today, ".",VARIANT,".regional_assoc.eps"), plot = last_plot())
  
  rm(temp, p1,
     temp_f, temp_f_ld,
     tempCHR, tempSTART, tempEND,
     VARIANT, tempVARIANTnr)
  
}

Getting data for rs10899970.

Subset required data.

Formatting association data.

Formating association data...
Processing -log10(p-values)...
Preparing association data...

Getting LD data.

All inputs are go!
Reading in association data...
Populations selected.
Calculating LD using rs10899970...

[1] "https://ldlink.nci.nih.gov/LDlinkRest/ldproxy?var=rs10899970&pop=EUR&r2_d=r2&token=c0f613f149ab"

 [9%] Downloaded 15979 bytes...
 [20%] Downloaded 32209 bytes...
 [30%] Downloaded 48575 bytes...
 [50%] Downloaded 81343 bytes...
 [60%] Downloaded 97745 bytes...
 [71%] Downloaded 114111 bytes...
 [81%] Downloaded 130513 bytes...
 [100%] Downloaded 160476 bytes...
Merging input association data with LD...

Plotting region surrounding rs10899970 on 10:44015716-45334720.

Plotting by...
snp rs10899970
Reading in association data
Generating Plot

Saving image for rs10899970.

Saving 7.29 x 4.51 in image

Additional regional plots

We want to create some regional association plots to combine with teh UCSC browser tracks, thus we need the exact same regions.

library(openxlsx)
add_list <- read.xlsx(paste0(TARGET_loc, "/Variants.xlsx"), sheet = "AdditionalPlots")

DT::datatable(add_list)

NA

library(RACER)
# Make directory for plots
ifelse(!dir.exists(file.path(PROJECT_loc, "/RACER")), 
       dir.create(file.path(PROJECT_loc, "/RACER")), 
       FALSE)

[1] FALSE

RACER_loc = paste0(PROJECT_loc,"/RACER")

variants_of_interest <- c(add_list$rsID)


for(VARIANT in variants_of_interest){
  cat(paste0("Getting data for ", VARIANT,".\n"))

  tempCHR <- subset(add_list, rsID == VARIANT)[,4]
  tempSTART <- subset(add_list, rsID == VARIANT)[,5]
  tempEND <- subset(add_list, rsID == VARIANT)[,6]
  tempNAME <- subset(add_list, rsID == VARIANT)[,3]

  cat("\nSubset required data.\n")
  temp <- subset(gwas_sumstats_racer, Chr == tempCHR & (Position >= tempSTART & Position <= tempEND))
  
  cat("\nFormatting association data.\n")
  temp_f = RACER::formatRACER(assoc_data = temp, chr_col = 3, pos_col = 4, p_col = 5)

  cat("\nGetting LD data.\n")
  temp_f_ld = RACER::ldRACER(assoc_data = temp_f, rs_col = 2, pops = "EUR", lead_snp = VARIANT)
  
  cat(paste0("\nPlotting region surrounding ", VARIANT," on ",tempCHR,":",tempSTART,"-",tempEND,".\n"))
  # source(paste0(PROJECT_loc, "/scripts/functions.R"))
  p1 <- singlePlotRACER2(assoc_data = temp_f_ld, 
                         chr = tempCHR, build = "hg19", 
                         plotby = "coord", snp_plot = VARIANT,
                         start_plot = tempSTART, end_plot = tempEND,
                         label_lead = FALSE, 
                         grey_colors = TRUE, gene_track_h = 3, gene_name_s = 1.75)
  
  print(p1)
  
  cat(paste0("Saving image for ", VARIANT,".\n"))
  ggsave(filename = paste0(RACER_loc, "/", tempNAME, ".", Today, ".",VARIANT,".",tempSTART,".",tempEND,".regional_assoc.png"), plot = last_plot())
  ggsave(filename = paste0(RACER_loc, "/", tempNAME, ".", Today, ".",VARIANT,".",tempSTART,".",tempEND,".regional_assoc.pdf"), plot = last_plot())
  ggsave(filename = paste0(RACER_loc, "/", tempNAME, ".", Today, ".",VARIANT,".",tempSTART,".",tempEND,".regional_assoc.eps"), plot = last_plot())

  rm(temp, p1,
     temp_f, temp_f_ld,
     tempCHR, tempSTART, tempEND,
     VARIANT, tempNAME)
  
}

Getting data for rs9633535.

Subset required data.

Formatting association data.

Formating association data...
Processing -log10(p-values)...
Preparing association data...

Getting LD data.

All inputs are go!
Reading in association data...
Populations selected.
Calculating LD using rs9633535...

[1] "https://ldlink.nci.nih.gov/LDlinkRest/ldproxy?var=rs9633535&pop=EUR&r2_d=r2&token=c0f613f149ab"

 [100%] Downloaded 49363 bytes...
Merging input association data with LD...

Plotting region surrounding rs9633535 on 10:63584853-63921073.

Plotting by...
coord
Reading in association data
Generating Plot

Saving image for rs9633535.

Saving 7.29 x 4.51 in image
Saving 7.29 x 4.51 in image
Saving 7.29 x 4.51 in image

Getting data for rs2854746.

Subset required data.

Formatting association data.

Formating association data...
Processing -log10(p-values)...
Preparing association data...

Getting LD data.

All inputs are go!
Reading in association data...
Populations selected.
Calculating LD using rs2854746...

[1] "https://ldlink.nci.nih.gov/LDlinkRest/ldproxy?var=rs2854746&pop=EUR&r2_d=r2&token=c0f613f149ab"

 [65%] Downloaded 71861 bytes...
 [80%] Downloaded 88263 bytes...
 [100%] Downloaded 109298 bytes...
Merging input association data with LD...

Plotting region surrounding rs2854746 on 7:45894617-46054070.

Plotting by...
coord
Reading in association data
Generating Plot

Saving image for rs2854746.

Saving 7.29 x 4.51 in image
Saving 7.29 x 4.51 in image
Saving 7.29 x 4.51 in image

Getting data for rs3844006.

Subset required data.

Formatting association data.

Formating association data...
Processing -log10(p-values)...
Preparing association data...

Getting LD data.

All inputs are go!
Reading in association data...
Populations selected.
Calculating LD using rs3844006...

[1] "https://ldlink.nci.nih.gov/LDlinkRest/ldproxy?var=rs3844006&pop=EUR&r2_d=r2&token=c0f613f149ab"

 [100%] Downloaded 31016 bytes...
Merging input association data with LD...

Plotting region surrounding rs3844006 on 6:131937915-132289374.

Plotting by...
coord
Reading in association data
Generating Plot

Saving image for rs3844006.

Saving 7.29 x 4.51 in image
Saving 7.29 x 4.51 in image
Saving 7.29 x 4.51 in image

Session information

---

Version:      v1.2.0
Last update:  2022-01-28
Written by:   Sander W. van der Laan (s.w.vanderlaan-2[at]umcutrecht.nl).
Description:  Script to create plot regional association plots.
Minimum requirements: R version 3.4.3 (2017-06-30) -- 'Single Candle', Mac OS X El Capitan

Changes log
* v1.2.0 Added in aditional regions.
* v1.1.0 Created PNG and PDF of top loci regions.
* v1.0.0 Initial version. 

---

sessionInfo()

R version 4.1.2 (2021-11-01)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Monterey 12.2

Matrix products: default
LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] grid      tools     stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] RACER_1.0.0          openxlsx_4.2.5       UpSetR_1.4.0         ggrepel_0.9.1        ggpubr_0.4.0         forestplot_2.0.1     checkmate_2.0.0     
 [8] magrittr_2.0.1       pheatmap_1.0.12      devtools_2.4.3       usethis_2.1.5        BlandAltmanLeh_0.3.1 tableone_0.13.0      haven_2.4.3         
[15] eeptools_1.2.4       DT_0.20              knitr_1.37           forcats_0.5.1        stringr_1.4.0        purrr_0.3.4          tibble_3.1.6        
[22] ggplot2_3.3.5        tidyverse_1.3.1      data.table_1.14.2    naniar_0.6.1         tidyr_1.1.4          dplyr_1.0.7          optparse_1.7.1      
[29] readr_2.1.1         

loaded via a namespace (and not attached):
  [1] rappdirs_0.3.3              rtracklayer_1.54.0          scattermore_0.7             coda_0.19-4                 ragg_1.2.1                 
  [6] SeuratObject_4.0.4          visdat_0.5.3                bit64_4.0.5                 irlba_2.3.5                 DelayedArray_0.20.0        
 [11] rpart_4.1-15                KEGGREST_1.34.0             RCurl_1.98-1.5              AnnotationFilter_1.18.0     generics_0.1.1             
 [16] BiocGenerics_0.40.0         GenomicFeatures_1.46.3      callr_3.7.0                 cowplot_1.1.1               RSQLite_2.2.9              
 [21] RANN_2.6.1                  future_1.23.0               bit_4.0.4                   tzdb_0.2.0                  spatstat.data_2.1-2        
 [26] xml2_1.3.3                  lubridate_1.8.0             httpuv_1.6.5                SummarizedExperiment_1.24.0 assertthat_0.2.1           
 [31] xfun_0.29                   jquerylib_0.1.4             hms_1.1.1                   evaluate_0.14               promises_1.2.0.1           
 [36] fansi_1.0.0                 restfulr_0.0.13             progress_1.2.2              dbplyr_2.1.1                readxl_1.3.1               
 [41] igraph_1.2.11               DBI_1.1.2                   htmlwidgets_1.5.4           spatstat.geom_2.3-1         stats4_4.1.2               
 [46] ellipsis_0.3.2              crosstalk_1.2.0             backports_1.4.1             survey_4.1-1                biomaRt_2.50.1             
 [51] deldir_1.0-6                MatrixGenerics_1.6.0        vctrs_0.3.8                 Biobase_2.54.0              remotes_2.4.2              
 [56] ensembldb_2.18.2            ROCR_1.0-11                 abind_1.4-5                 cachem_1.0.6                withr_2.4.3                
 [61] vcd_1.4-9                   sctransform_0.3.2           GenomicAlignments_1.30.0    prettyunits_1.1.1           getopt_1.20.3              
 [66] goftest_1.2-3               cluster_2.1.2               lazyeval_0.2.2              crayon_1.4.2                labeling_0.4.2             
 [71] pkgconfig_2.0.3             GenomeInfoDb_1.30.0         nlme_3.1-153                pkgload_1.2.4               ProtGenerics_1.26.0        
 [76] rlang_0.4.12                globals_0.14.0              lifecycle_1.0.1             miniUI_0.1.1.1              filelock_1.0.2             
 [81] BiocFileCache_2.2.0         modelr_0.1.8                cellranger_1.1.0            rprojroot_2.0.2             polyclip_1.10-0            
 [86] matrixStats_0.61.0          lmtest_0.9-39               Matrix_1.4-0                carData_3.0-5               boot_1.3-28                
 [91] zoo_1.8-9                   reprex_2.0.1                ggridges_0.5.3              processx_3.5.2              png_0.1-7                  
 [96] viridisLite_0.4.0           rjson_0.2.21                bitops_1.0-7                KernSmooth_2.23-20          pander_0.6.4               
[101] Biostrings_2.62.0           blob_1.2.2                  maptools_1.1-2              arm_1.12-2                  parallelly_1.30.0          
[106] rstatix_0.7.0               S4Vectors_0.32.3            ggsignif_0.6.3              scales_1.1.1                memoise_2.0.1              
[111] plyr_1.8.6                  ica_1.0-2                   zlibbioc_1.40.0             compiler_4.1.2              BiocIO_1.4.0               
[116] RColorBrewer_1.1-2          lme4_1.1-27.1               fitdistrplus_1.1-6          Rsamtools_2.10.0            cli_3.1.0                  
[121] XVector_0.34.0              listenv_0.8.0               patchwork_1.1.0.9000        pbapply_1.5-0               ps_1.6.0                   
[126] MASS_7.3-54                 mgcv_1.8-38                 tidyselect_1.1.1            stringi_1.7.6               textshaping_0.3.6          
[131] mitools_2.4                 yaml_2.2.1                  sass_0.4.0                  future.apply_1.8.1          parallel_4.1.2             
[136] rstudioapi_0.13             foreign_0.8-81              gridExtra_2.3               farver_2.1.0                Rtsne_0.15                 
[141] digest_0.6.29               BiocManager_1.30.16         shiny_1.7.1                 Rcpp_1.0.7                  GenomicRanges_1.46.1       
[146] car_3.0-12                  broom_0.7.11                later_1.3.0                 RcppAnnoy_0.0.19            httr_1.4.2                 
[151] AnnotationDbi_1.56.2        colorspace_2.0-2            rvest_1.0.2                 XML_3.99-0.8                fs_1.5.2                   
[156] tensor_1.5                  reticulate_1.22             IRanges_2.28.0              splines_4.1.2               uwot_0.1.11                
[161] spatstat.utils_2.3-0        sp_1.4-6                    systemfonts_1.0.3           plotly_4.10.0               sessioninfo_1.2.2          
[166] xtable_1.8-4                nloptr_1.2.2.3              jsonlite_1.7.2              testthat_3.1.1              R6_2.5.1                   
[171] pillar_1.6.4                htmltools_0.5.2             mime_0.12                   minqa_1.2.4                 glue_1.6.0                 
[176] fastmap_1.1.0               BiocParallel_1.28.3         codetools_0.2-18            pkgbuild_1.3.1              utf8_1.2.2                 
[181] bslib_0.3.1                 lattice_0.20-45             spatstat.sparse_2.1-0       curl_4.3.2                  leiden_0.3.9               
[186] zip_2.2.0                   survival_3.2-13             rmarkdown_2.11              desc_1.4.0                  munsell_0.5.0              
[191] GenomeInfoDbData_1.2.7      reshape2_1.4.4              gtable_0.3.0                spatstat.core_2.3-2        

Saving environment

save.image(paste0(PROJECT_loc, "/",Today,".",PROJECTNAME,".RegionalAssociationPlots.RData"))

© 1979-2022 Sander W. van der Laan | s.w.vanderlaan[at]gmail.com | swvanderlaan.github.io.